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We extend the recently introduced phaseless auxiliary-field quantum Monte Carlo (QMC) ap- 
proach to any single-particle basis, and apply it to molecular systems with Gaussian basis sets. 
QMC methods in general scale favorably with system size, as a low power. A QMC approach with 
auxiliary fields in principle allows an exact solution of the Schrodinger equation in the chosen basis. 
However, the well-known sign/phase problem causes the statistical noise to increase exponentially. 
The phaseless method controls this problem by constraining the paths in the auxiliary-field path 
integrals with an approximate phase condition that depends on a trial wave function. In the present 
calculations, the trial wave function is a single Slater determinant from a Hartree-Fock calculation. 
The calculated all-electron total energies show typical systematic errors of no more than a few 
milli-Hartrees compared to exact results. At equilibrium geometries in the molecules we studied, 
this accuracy is roughly comparable to that of coupled-cluster with single and double excitations 
and with non- iterative triples, CCSD(T). For stretched bonds in H2O, our method exhibits better 
overall accuracy and a more uniform behavior than CCSD(T). 

PACS numbers: 02.70.Ss, 71.15.-m, 31.25.-v 
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I. INTRODUCTION 

Obtaining the solution of the Schrodinger equation, 
even within a finite Hilbert space, is an important goal 
in quantum chemistry and in condensed matter physics, 
both in its own right and for calibrating approximate 
methods. In quantum chemistry, a variety of approaches 
have been developed over the last fifty years, which range 
in quality from the mean-field Hartree-Fock (HF) to the 
exact full configuration interaction (FCI) solution. In 
between these two limits, a hierarchy of approximations 
have been developed which improve upon the Hartree- 
Fock solution at the expense of increasing computational 
cost [1]. 

The FCI method is quite formidable computationally, 
because its cost grows exponentially with the number of 
electrons and basis functions. Recently, the density ma- 
trix renormahzation group (DMRG) approach was intro- 
duced as a method to potentially obtain near-FCI quality 
solutions of the Schrodinger equation for quantum chem- 
ical Hamiltonians [2-4]. The results so far are promising. 
For example, DMRG was applied recently [3] to study 
the water molecule using a double-zeta atomic-natural- 
orbital (ANO) basis (41 basis functions), which is cur- 
rently intractable with FCI. 

Among the approximate approaches, coupled cluster 
(CC) methods have been the most successful, especially 
at equilibrium geometries [5, 6]. However, standard 
CC methods are non- variational and their computational 
cost grows rather steeply with the system size as well. For 
example, the most popular among them, CC with sin- 
gle and double excitations and with non-iterative triples 
[CCSD(T)], has an N"^ scaling with the size of the basis. 

Quantum Monte Carlo (QMC) methods are an attrac- 
tive means for a non-perturbative and explicit treatment 
of the interacting many-fermion system. These methods 



tend to have favorable scaling with system size, often as a 
low power. The most established QMC method, the real 
space fixed-node diffusion Monte Carlo (DMC) method, 
has been applied successfully to calculate many proper- 
ties of solids [7] and molecules [8] . 

Recently, an alternative and complementary QMC 
method has been introduced [9] for realistic electronic 
Hamiltonians, based on a field-theoretic approach with 
auxiliary fields. The central idea of auxiliary-field 
(AF) QMC methods is the mapping of the interacting 
many-body problem into a linear combination of non- 
interacting problems coupled to external AFs. Averaging 
over different AF configurations is performed by Monte 
Carlo techniques. The basic formalism of AF QMC meth- 
ods has mostly been applied to lattice models of strongly 
interacting systems [10, 11]. In these models, the sim- 
plified form of the two-body interactions makes possible 
an efficient mapping with real AFs. As a result, "only" 
a sign problem [12] occurs. Constrained path methods 
[12, 13] have been developed to approximately control 
the sign problem, which have been shown to be quite 
accurate. 

The potential of AF QMC methods for realistic Hamil- 
tonians has long been recognized and pursued [14-16]. 
However, wide and general applications were not real- 
ized because of a lack of control for the phase problem. 
Except for special cases, the two-body electronic interac- 
tions lead to complex AFs. The many-dimensional inte- 
gration over complex variables in turn leads to an expo- 
nentially increasing noise and therefore breakdown of the 
method. The phaseless AF QMC method [9] was intro- 
duced to control the phase problem in an approximate 
manner. 

In Ref. [9] , a general framework was developed on how 
to use importance sampling of the determinants to deal 
with the complex phase. An approximate constraint was 
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then formulated with a trial wave function to constrain 
the paths in AF space and eliminate the growth of noise. 
Because of the constraint, the calculated ground state 
energy is no longer exact and is not variational. In ap- 
plications to several sp-bonded materials [9, 17, 18], it 
was shown that the method, with a plane-wave basis and 
simple trial wave functions, gave accurate results for sys- 
tems from simple atoms to large supercells. The phase- 
less AF QMC method was also recently applied to the 
strongly correlated transition metal oxide molecules TiO 
and MnO [19], again yielding results in good agreement 
with experiment using simple mean-field trial wave func- 
tions. 

The success of the phaseless AF approach in solid-state 
applications with plane- waves has motivated us to extend 
it to a generic one-particle basis, targeting in particu- 
lar quantum chemistry problems. Applying the new AF 
QMC method to such problems is very appealing, espe- 
cially with the integral role basis sets play in quantum 
chemistry methods. This method allows QMC calcula- 
tions using any choice of one-particle basis. It provides a 
framework for general many-body calculations which can 
import straightforwardly many of the well-established 
technologies from more traditional approaches. By 
importance-sampled random walks, the phaseless AF 
QMC method obtains the many-body ground state with 
an ensemble of loosely coupled independent-particle cal- 
culations which propagate with imaginary-time. The 
ground-state wave function is represented by a linear 
combination of non-orthogonal Slater determinants. 

In this paper we report the first results from this ef- 
fort. We discuss aspects of the new method when im- 
plemented with Gaussian basis sets. We present bench- 
mark results on sp-bonded atoms and molecules. All- 
electron total energies are calculated for various first- 
row atoms and molecules, and compared with FCI and 
DMRG results where available or with CCSD(T) other- 
wise. The behavior of the method is also studied as a 
function of basis size, including an ANO basis [20, 21] 
in H2O and O2, each with 92 basis fmictions. Wo then 
test the method away from equilibrium geometries, cal- 
culating the equilibrium bond length and also stretching 
the bonds in the water molecule by a factor between 1 
and 8 within a cc-pVDZ basis [22]. In all our calculations 
the trial wave function is taken to be a single Slater de- 
terminant from Hartree-Fock, either restricted (RHF) or 
unrestricted (UHF). Their effect on the accuracy is exam- 
ined. The calculated energies with the optimal HF wave 
function typically show systematic errors of no more than 
a few milli-Hartree (mEh) compared to exact results. For 
stretched bonds in H2O, our method exhibits better over- 
all accuracy and more uniform behavior than CCSD(T). 

An additional benefit of the present AF QMC imple- 
mentation is its value for algorithmic studies. Compared 
to the plane-wave algorithm, the generic basis AF QMC 
method typically has much more effective basis sets. As 
a result, the basis size is often much smaller and the 
method is significantly cheaper computationally. Fur- 



ther, direct comparisons can be more easily made be- 
tween the QMC results and those obtained with more 
established correlated methods, as indicated above. Ex- 
tra ingredients such as pseudopotentials can be removed, 
with the only remaining uncertainty being the systematic 
error from the phaseless approximation. In this study, we 
take advantage of this feature to make detailed bench- 
mark studies of the systematic error. 

The rest of the paper is organized as follows. The 
phaseless AF QMC method is first briefly reviewed in 
the next section. We then describe its implementation 
with the Gaussian basis in Sec. Ill, together with some 
results to illustrate the behavior and characteristics of 
this approach. Section IV gives some of the details of 
the computational procedures. In Sec. V and Sec. VI, we 
present the results of our calculations, including calcula- 
tions of the II2O molecule both at the equilibrium geom- 
etry and with stretched bond lengths, and total energies 
and energy differences (binding, ionization, and electron 
affinity) for various other first-row atoms and diatomic 
molecules. Finally in Sec. VII we make concluding re- 
marks, together with a brief discussion of possibilities for 
further improvement of the new algorithm. 

II. FORMALISM 

The full electronic Hamiltonian for a many-fermion 

system with two-body interactions can be written in any 
orthonormal one-particle basis in the general form 

N I ^ 

H = H,+H2=J2 TiAcj + 2 E VijkicWjCkCi, (1) 

i,j i,j,k,l 

where N is the size of the chosen one-particle basis, and 
c| and c, are the corresponding creation and annihilation 
operators. Both the one-body (Tij) and two-body (Vijki) 
matrix elements are known. 

The auxiliary-field quantum Monte Carlo method ob- 
tains the ground state \^g) of H, by projecting from a 
trial wave function \^t) using the imaginary- time prop- 
agator e~^^: 

\-^g) oc lim (e-^-^)" |*t) • (2) 

n— >oo 

The trial wave function \^t), which should have a non 
zero overlap with the exact ground state, is assumed to 
be in the form of a single determinant or a linear combi- 
nation of Slater determinants. 

The timcstcp, r, is chosen to be small enough so that 
Hi and H2 in the propagator can be accurately separated 
with the Trotter decomposition: 

The action of the propagator e~'^^^ on a Slater deter- 
minant yields another determinant. This is not the case 
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with e^'^ , which however can be written as an inte- 
gral of one-body operators using a Hubbard-Stratonovich 
(HS) transformation [23]: 




where the one-body operators Va can be defined gener- 
ally for any two-body operator by writing the latter in 
a quadratic form, such as H2 = "^X^aCa'^'ai with C,a 
a real number. Monte Carlo methods are very efficient 
at evaluating multi-dimensional integrals as in Eq. (4). 
For example, the projection of Eq. (2) can be realized 
iteratively: An ensemble of Slater determinants { |0) } 
are initialized to the trial wave function |^'t), which are 
then propagated to a new ensemble { |(/)') } using e""^^ 
of Eq. (3), and so on, until convergence is reached. For 
each Slater determinant in the ensemble, the two-body 
part in Eq. (4) can be applied stochastically by sampling 
an AF configuration, {era}. In the standard AF QMC 
approach, the projection is often done [12] as a path in- 
tegral with paths of fixed length in imaginary time, using 
a Metropolis or heat-bath algorithm. 

This straightforward approach, however, will generally 
lead to an exponential increase in the statistical fluctu- 
ations with the projection time. The source of the fluc- 
tuations is that the one-body operators v = {va} are 
generally complex, since usually Ca cannot all be made 
positive in Eq. (4). As a result, the orbitals in will 
become complex as the projection proceeds. For large 
projection time, the phase of each \<f>) becomes random, 
and the MC representation of |4'g) becomes dominated 
by noise. This leads to the phase problem and the diver- 
gence of the fluctuations. The phase problem is of the 
same origin as the sign problem that occurs when the one- 
body operators v are real, but is more severe because for 
each |(^), instead of a and —\4>) symmetry [13], there 
is now an infinite set {e*^|^),^ e [0, 27r)} among which 
the Monte Carlo sampling cannot distinguish. 

We used the phaseless auxiliary-field QMC method to 
control the phase problem. In order to implement a 
phaseless constraint, the method recasts the imaginary- 
time path integral as a branching random walk in Slater- 
determinant space [9, 13]. It uses a trial wave function 
I^t) and a comp/ea; importance function, (^Tltp), to con- 
struct phaseless random walkers, / {^t\4') , which are 
invariant under a phase gauge transformation. The re- 
sulting two-dimensional diffusion process in the complex 
plane of the overlap (^I/tI^) is then approximated as a 
diffusion process in one dimension. 

As mentioned earlier, the ground-state energy com- 
puted with the so-called mixed estimate is approximate 
and not variational in the phaseless method. The error 
depends on |^'t), vanishing when |5't) is exact. This 
is the only error in the method that cannot be elimi- 
nated systematically. The method has been applied to 
atoms, molecules, and simple solids, using a plane-wave 
basis and pseudopotentials, and has proved very success- 



ful [9, 17, 19]. 

III. IMPLEMENTATION USING A LOCALIZED 
BASIS 

With Gaussian basis sets, the matrix elements of inter- 
est to QMC and the overlap matrix can be imported from 
quantum chemistry programs. It is more convenient to 
use orthonormal basis functions, which ensure the usual 
commutation relations of the creation and destruction 
operators. To set up the Hamiltonian in Eq. (1), we first 
transform the non-orthogonal basis into an orthogonal- 
ized basis set. The one- and two-body matrix elements, 
Tij and Vijki , are then expressed with respect to this ba- 
sis via the transformation, based on the original matrix 
elements. 

To carry out the HS transformation of Eq. (4), we 
first map the two-body interaction matrix Vijki into a 
real, symmetric supermatrix Vf^[i^i]^^[j^k] where fijV = 
1, . . . , iV^. This is then expressed in terms of its eigenval- 
ues (— Aa) and eigenvectors X^: V^^i, = — J2a ^aX° X^. 

The two-body operator H2 of Eq. (1) can be written 
as the sum of a one-body operator H[ and a two-body 
operator H2 ■ The latter can be further expressed in terms 
of the eigenvectors of V^^i, as 

a 

where the one-body operator Va is, 

i,l 

Note that — Va for real basis functions, as is the case 
here. The form in Eq. (5) is now amenable to the HS 
transformation of Eq. (4). 

In this case the number of the HS fields is equal to 
the number of non-zero eigenvalues (— Aq) of the sym- 
metrized two-body supermatrix. Other forms of decom- 
position and HS transformation are possible, and their ef- 
ficiency and effectiveness with a constraint can vary. (For 
example, using a modified Cholesky decomposition will 
eliminate the need to diagonalize the supermatrix and 
potentially reduces the computational cost of synthesiz- 
ing the matrix. However, this will generally lead to a 
larger number of HS fields, and hence increase the QMC 
computational cost.) We will defer further analysis of 
this freedom to a future publication. 

To illustrate the method, we use the H2O molecule in a 
minimal ST0-6G basis. This problem is small enough to 
permit detailed comparison with exact diagonalization. 
Starting from HF solution, the projection in Eq. (2) was 

applied using: i) the exact many-body propagator e~^^ 
acting on the wave function expanded in terms of the 

exact cigcnstatcs of the many-body Hamiltonian H, ii) 
the free projection QMC method (no constraint imposed) 
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FIG. 1: The projection of the ground state of H2O/STO-6G. 
Starting from the HF trial wave function, the total energy is 
plotted vs. the projection time /3 = nr, as obtained using 
the free projection QMC (exact but suffers from the phase 
problem) and the phaseless approximation. For comparison, 
we show also the exact projection which is obtained form di- 
agonalization of the many-body Hamiltonian. The phaseless 
QMC method tracks smoothly the exact projection with no 
sign of the uncontrolled fluctuations in free projection. 



which is exact but suffers from the phase problem, and 
iii) the phaseless QMC which is approximate but controls 
the phase problem. Figure 1 plots the results of total- 
energy mixed-estimator [9] , corresponding to these three 
calculations. Up to a projection time of /3 ^ 2 E 
both QMC methods show essentially the same conver 
gence of the total-energy. For large projection times, the 
free-projection starts showing the phase problem in the 
form of large fluctuations. Around /3 « 23 E^^ these 
fluctuations diverge on the scale of the plot. For larger 
number of particles or larger basis sets, this divergence 
would have occurred at much earlier projection times (3. 
The phaseless QMC method, by contrast, follows the ex- 
act projection with finite variance as the random walk 
proceeds. 

In our implementation, we found it advantageous to 
first subtract the mean-field "background" from Va be- 
fore applying the Hubbard-Stratonovich transformation. 
In this case, Eq. (5) is replaced by. 



1 

h ' 



(7) 



where (va) — (^tIwoI^'t) and H'l is the one-body oper- 
ator: 

H'^ = -\ Xl(Wa(Wa) +Va(Wa) " (8) 



Equations (5) and (7) are equivalent. When the approx- 
imate phaseless projection is imposed, however, using 
Eq. (7) before applying the HS transformation leads to 
an improved behavior. 
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FIG. 2: QMC Trotter time-step dependence. The total en- 
ergy of the H2O/STO-6G molecule from phaseless QMC is 
shown vs. the time-step size r, using the HS transforma- 
tions of Eq. (5) (labeled as QMC/NS) and of Eq. (7) (labeled 
as QMC/WS, in which the mean-field "background" is sub- 
tracted). Results from free-projection (exact but suffers from 
the phase problem) are also shown, together with that of the 
exact FCI. Points are connected by straight lines for clarity. 
The insert shows a closer look at small values of r, together 
with linear fits of the QMC data. 



We illustrate this point in Fig. 2. Using the same STO- 
6G minimal basis, we plot the Trotter time-step conver- 
gence of the total energy of the H2O molecule, calculated 
with the HS transformations of Eq. (5) and Eq. (7), for 
both the phaseless QMC method and for (exact) free- 
projection. For comparison, the ground state energy ob- 
tained from exact diagonalization is also shown. The 
Trotter extrapolated QMC energy (see inset of Fig. 2) is 
-75.72725(4) E,, with the HS decomposition of Eq. (5) 
and -75.72849(5) E/^ with Eq. (7). Monte Carlo sta- 
tistical error bars are on the last digit and are shown 
in parentheses. For comparison, the exact ground-state 
energy is -75.72799 E,,. 

Similar results and analysis have also been presented 
in phaseless AF QMC calculations in boson systems [24] . 
We comment that the issue in Fig. 2 is different from 
the shifted contour approach [15], dispite formal similar- 
ities. With the importance sampling scheme of Ref. [9], 
Eqs. (5) and (7) would both be exact and give the same 
results under free projection, since the mean-field shift is 
automatically applied via the importance sampling trans- 
formation regardless of which form of H'^ is used [25]. 
The different behaviors discussed above arise only be- 
cause of the imposition of the phase projection to one- 
dimension [9] , in which the substraction of the mean-field 
background helps to reduce the "rotation" of the random 
walkers in the complex (5'T|(/')-plane and thus the sever- 
ity of the projection. 
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IV. COMPUTATIONAL DETAILS 

We apply our method to several atomic and molecular 
systems, which were chosen primarily for benchmarking 
purposes and which have all been previously studied us- 
ing well-established correlated methods such as coupled 
cluster (CC), full configuration interaction (FCI), or den- 
sity matrix renormalization group (DMRG) methods. In 
these and in our QMC calculations, the core and valence 
electrons are fully correlated. 

All the QMC calculations reported below used a single 
Hartrcc-Fock Slater determinant as the trial wave func- 
tion. For all systems, QMC was done with the best 
variational single determinant, namely the unrestricted 
Hartrec-Fock (UHF) solution. In addition, we have also 
tested the restricted Hartree-Fock (RHF) solution as the 
trial wave function in some cases. No additional opti- 
mization or tuning of the trial wave function was per- 
formed beyond the mean-field Hartree-Fock calculation. 

The QMC code is interfaced with GAMESS[26] and 
NWChem[27] to import the Gaussian one-electron and 
two-electron matrix elements, the overlap matrix, and 
the trial wave function. All of our calculations are done 
using the spherical harmonics representation of basis 
functions. 

We performed the FCI calculations using MOL- 
PRO [28, 29] and the coupled cluster calcuations using 
NWChem [27] . Unless otherwise noted, the coupled clus- 
ter calculations for atoms and molecules with a singlet 
state are of the type RIIF-RCCSD(T) i.e. based on the 
RHF reference state, while those for a multiplet state are 
of the unrestricted type, UHF-UCCSD(T) based on the 
UHF solution. 



V. EQUILIBRIUM GEOMETRIES 

As a first test of our method, we calculate in this sec- 
tion total energies and binding energies for some well- 
studied systems at their equilibrium geometries. The 
first subsection contains results for the water molecule 
in the vicinity of the equilibrium structure and the O2 
dimer at equilibrium bond length. The second subsection 
contains results for other first-row atoms and diatomic 
molecules, including total ground-state energies, binding 
energy, convergence with basis sets, atomic ionization po- 
tentials, and electron affinities. In the next section, we 
test the method by studying bond stretching of the water 
molecule. 



A. H2O and O2 

The water molecule has a long benchmarking history 
[3, 30-33]. Table I presents results calculated by HF, 
CCSD(T), FCI, and the present QMC method. We show 
the all-electron total energies of the individual atoms 
H and O, and the total and binding energies of H2O. 



TABLE I: The binding energy (BE) of H2O calculated with 
four basis sets: STO-6G, cc-pVDZ, and double- and triple- 
zeta ANO [21]. Also shown are all-electron total energies 
for O, H, and H2O. Monto Carlo statistical errors arc in the 
last digit and arc shown in parentheses. All energies are in 
Hartrees. The FCI value for H20/cc-pVDZ is from Ref. [32], 
while the DMRG energy of H20 and the FCI energy of O 
within the double-zeta ANO are from Ref. [3] . 



UHF CCSD(T) FCI/DMRG QMC 

ST0-6G 

H -0.471039 -0.471039 -0.471039 -0.471039 

O -74.516 816 -74.516 816 -74.516 816 -74.516 816 

H2O -75.676 506 -75.727 931 -75.727 991 -75.728 5(1) 

BE 0.217612 0.269 037 0.269 097 0.269 5(1) 

cc-pVDZ 

H -0.499 278 -0.499 278 -0.499 278 -0.499 278 

O -74.792166 -74.911552 -74.911744 -74.909 6(1) 

H2O -76.024 039 -76.241 201 -76.241 860 -76.242 4(2) 

BE 0.233 317 0.331093 0.331560 0.334 3(2) 

DZ ANO 

H -0.499 944 -0.499 944 -0.499 944 -0.499 944 

O -74.816 273 -74.961956 -74.962 350 -74.959 6(1) 

H2O -76.057 621 -76.314141 -76.314 715 -76.316 3(6) 

BE 0.241460 0.352 959 0.352 477 0.356 8(7) 

TZ ANO 

H -0.499 973 -0.499 973 -0.499 973 

O -74.818 648 -75.000129 -74.9970(4) 

H2O -76.060 589 -76.367 528 -76.370(1) 

BE 0.241995 0.367 453 0.373(1) 



Four different basis sets were used: minimal ST0-6G, cc- 
pVDZ [22] , and the double zeta (DZ) and triple-zeta (TZ) 
ANO bases of Widmark, Malmqvist, and Roos [20, 21]. 
For H2O, these have 7, 24, 41, and 92 basis functions, re- 
spectively. The H2O geometries were as follows: results 
using the minimal ST0-6G correspond to positions (in 
Bohr) of 0(0,0,0) and H(0, ±1.515263, -1.058898) with 
C2v symmetry, while the cc-pVDZ results used the ge- 
ometry of Ref. [32]. The DZ ANO and TZ ANO results 
were obtained with the geometry given in Ref. [3]. 

The QMC calculations used UHF trial wave functions. 
(For H2O, the RHF and UHF solutions are identical at 
the equilibrium bond length.) The FCI total energies of 
H2O/CC-PVDZ arc those of Ref. [32]. FCI resuhs for H2O 
are not within reach for the DZ ANO basis, but DMRG 
results are available from a recent study [3], which are 
shown under FCI. The FCI energy of 0/DZ ANO is also 
from Ref. [3]. Neither FCI nor DMRG is within reach 
presently for the TZ ANO basis. 

We see that the agreement between the QMC, 
CCSD(T), and FCI results is quite good. As mentioned, 
the calculated groimd-statc energies in the present phase- 
less QMC method are not variational, which can be seen 
in the results in H2O, for example. QMC binding energy 
results of H2O overestimate the exact values by 0.4(1), 
2.7(2), 4.4(6) mEh, for the ST0-6G, cc-pVDZ,' and DZ 
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FIG. 3: Calculation of the equilibrium bond length of H2O 
within the cc-pVDZ basis and DZ ANO basis, using CCSD(T) 
and the present QMC method. All-electron total energies (in 
Hartrees) are shown as a function of the O-H bond length 
(in units of the experimental equilibrium value, Re), with the 
angle between the two O-H bonds fixed at the experimental 
value. Monte Carlo statistical error bars are indicated. The 
lines are based on a polynomial fit to the data. 



TABLE II: Binding energy (BE) of O2 using cc-pVDZ, 
double- and triple-zeta ANO [21] basis sets. The calculated 
all-electron total energies are also shown. The experimental 
bond length of Re = 1.2074 A was used. O atom results are 
in Table I. Monte Carlo statistical errors are in the last digit 
and are shown in parentheses. All energies are in Hartrees. 



UHF UCCSD(T) 



QMC 



cc-pVDZ 
O2 
BE 



-149.627 780 
0.043 448 



DZ ANO 

O2 

BE 

TZ ANO 

O2 

BE 



149.679 669 
0.047123 



149.989 319 
0.166 215 



150.095 668 
0.171 756 



-149.982 6(3) 
0.163 5(3) 



-150.090 7(6) 
0.1715(7) 



149.689 251 -150.187 965 -150.186(2) 
0.051955 0.187 706 0.192(2) 



ANO basis, respectively. The CCSD(T) method, which 
is highly accurate at this geometry, has errors of 0.06, 
0.47, 0.48 mEh, respectively. In the TZ ANO basis, 
CCSD(T) and QMC yield binding energies of 0.3675 and 
0.373(1) Eh, respectively. 

As mentioned, the DZ ANO and TZ ANO calcula- 
tions of H2O were for the geometry in Ref. [3]. Using 
CCSD(T), we verified that the experimental H2O geom- 
etry would lower the molecular total energy by 2.7 mEh 
and 4.0 mEh for the two basis sets. This would increase 
the binding energy by the same amount, and bring the 
CCSD(T) binding energy with the TZ ANO basis to good 
agreement with 0.370 Eh, the basis extrapolated value of 
Feller and Peterson [34]. For comparison, the experimen- 



tal binding energy of H2O is 0.3707 Eh (with the zero- 
point energy removed) [34]. 

Figure 3 shows a comparison of QMC and CCSD(T) re- 
sults in the vicinity of the H2O equilibrium bond length, 
using cc-pVDZ and ANO double zeta basis sets. The an- 
gle between the two O-H bonds is fixed at the experimen- 
tal angle of 104.4798 degrees, and we varied R/Rg where 
Re = 1.81 Bohr is the experimental bond length [35]. For 
the cc-pVDZ basis, the CCSD(T) and QMC equihbrium 
bond lengths calculated from these curves are 1.007 Re 
and 1.007(2) Re, respectively. The corresponding values 
using the ANO DZ basis set are 1.003 Re and 1.006(4) Re. 

We present results for the O2 molecule in Table II. The 
experimental equilibrium bond length of Re = 1.2074 A 
is used for all of these calculations. Results are shown 
for the three larger basis sets used in the above II2O cal- 
culations. These include the cc-pVDZ and the DZ and 
TZ ANO basis sets, which give 28, 46, and 92 basis func- 
tions for O2, respectively. QMC and CCSD(T) results 
are again in good agreement. In the TZ ANO basis, 
the CCSD(T) and QMC binding energies are 0.188 and 
0.192(2) Eh, respectively. Our QMC value is in excellent 
agreement with our previous binding energy 0.191(4) cal- 
culated using pseudopotentials and planewaves [19]. Also 
for comparison, the basis extrapolated binding energy 
at the CCSD(T) level is 0.191 Eh and the experimental 
value is 0.192 Eh (with the zero-point energy removed) 
[34]. 



B. First-row atoms and diatomic molecules 

1. All-electron total ground-state energies 

In Table III, we show the total energies of various first- 
row atoms and diatomic molecules at the equilibrium ge- 
ometry as calculated using RHF, UHF, CCSD(T), and 
also the FCI energies where available. For reference, we 
show also some available E'^^., the estimated experimen- 
tal non relativistic, infinite nuclear mass energy [36-38]. 
These represent the estimated exact results at the infinite 
basis size limit, and are thus significantly lower than the 
CCSD(T), FCI, or QMC energies because of the small 
basis size chosen in these benchmark calculations. 

For singlets where both RHF and UHF solutions exist, 
QMC calculations were done with each as a trial wave 
function, and both results are shown in Table III. In 
such cases, QMC with the best variational single Slater 
determinant as trial wave function, namely QMC/UHF, 
appears to always give better energy values. We discuss 
this effect further in Sec. VI with stretched bonds in H2O. 

Our QMC energies are generally within a few mEh of 
the FCI or the CCSD(T) energies. It is interesting to 
note the cases of the berylium atom and molecule. The 
Be atom, which has a near 2s-2p degeneracy, has often 
been used as a benchmark in Green's function or diffusion 
Monte Carlo (DMC) studies [40, 41]. Optimized Slater- 
Jastrow trial wave functions with a single determinant 
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TABLE III: The calculated all-electron total energies of various first-row atoms and diatomic molecules in their ground states. 
We show the RHF, UHF, CCSD(T) energies, and the FCI results where available. The basis set and the number of basis 
functions are shown in the third column. Our QMC results based on the RHF and UHF trial wave functions are shown in two 

columns. MC statistical errors arc in the last digit and shown in parentheses. For completeness, wc show in the last coluinn 
some of the available c;stiiiiate(l experimental, iion-relativistic:. infinite inicjlear mass energy (i.e.. the estimated exact I'esnlt at 
the iiifiiutt^ basis size limit). .Vll t'lit-rgies are in llartret's. 







uasisi ^ ) 


RHT7 


U LLC 










poo 


H2 


0.741 


cc-pVDZ(lO) 


-1.128 714 




-1.163 411 


-1.163 411 


-1.163 60(1) 




-1.1744 


Li 




cc-pVDZ(14) 




-7.432 421 


-7.432 638 


-7.432 638 




-7.432 633(3) 


-7.478 1 


Li2 


2.673 


cc-pVDZ(28) 


-14.869 499 


-14.870128 


-14.901386 


-14.901 392 


-14.902 5(1) 


-14.900 35(9) 


-14.995 4 


Be 




cc-pVDZ(14) 


-14.572 338 


-14.572 611 


-14.617407 


-14.617409 


-14.617 26(7) 


-14.617 63(8) 


-14.6674 


Bea 


2.45 


ec-pVDZ(28) 


-29.132 211 


-29.154 535 


-29.234 246 




-29.2313(2) 


-29.2341(2) 


-29.338 5 


Be 




cc-pVTZ(30) 


-14.572 874 


-14.573 183 


-14.623 790 


-14.623 810 


-14.622 4(2) 


-14.622 8(2) 


-14.667 4 


Bea 


2.45 


cc-pVTZ(60) 


-29.133 688 


-29.156 165 


-29.253 734 




-29.248 8(5) 


-29.251 1(3) 


-29.338 5 


N 




cc-pVDZ(14) 




-54.391 115 


-54.479 944 


-54.480 115 




-54.479 56(7) 


-54.589 3 


N2 


1.094 


cc-pVDZ(28) 


-108.954 553 




-109.278 722 




-109.2816(3) 




-109.542 3 


F 




cc-pVDZ(14) 




-99.375 240 


-99.529 322 


-99.529 518 




-99.5281(3) 


-99.734 


F2 


1.412 


cc-pVDZ(28) 


-198.685 670 


-198.695 746 


-199.101 151 




-199.100 3(5) 


-199.102 4(4) 




BH 


1.256 


cc-pVDZ(19) 


-25.125 188 


-25.131 227 


-25.215 917 


-25.216 401 


-25.217 5(3) 


-25.2151(2) 




CH+ 


1.146 


cc-pVDZ(19) 


-37.900 480 


-37.909 823 


-38.003 207 


-38.003 712" 


-38.006 9(3) 


-38.000 5(4) 




NH 


1.056 


cc-pVDZ(19) 




-54.966 091 


-55.093 298 


-55.093 721'' 




-55.093 18(7) 




OH+ 


1.032 


6-31G"(19) 




-74.973 345 


-75.093 371 






-75.093 08(9) 




HF 


0.920 


cc-pVDZ(19) 


-100.019 280 




-100.230 098 




-100.231 2(1) 







"Rcf. [39] 
''Ref. [39] 



are known to lead to significant errors, of ~ 10 mE/j, 

in the calculated fixed- node ground-state energy [40]. 
(This error is largely removed when multi-determinant 
trial wave functions are used to account for the near de- 
generacy. For example, DMC with an optimized Slater- 
Jastrow trial wave function using four determinants is 
extremely accurate [40].) In the present QMC method, 
the phasclcss approximation using a single Slater deter- 
minant appears to be significantly less dependent on the 
trial wave function, with systematic errors of ^ 1 mE^ or 
less. Of course, the DMC calculations work in real config- 
uration space and thus have the advantage of no finite- 
basis errors. Measured as a fraction of the correlation 
energy, the present QMC/UHF still has a substantially 
smaller systematic error, of about 2 Be/cc-pVTZ. 

Correspondingly, the Be2 dimer is a notorious case in 
quantum chemistry [42]. A DMC calculation using op- 
timized single-determinant trial wave functions did not 
obtain binding [43]. Our QMC/UHF total energies are 
within ~ 2 mEh of CCSD(T), and the resulting binding 
energy with the larger cc-pVTZ basis is 5.5(4) mE/i, in 
good agreement with the CCSD(T) value of 6.1 mE/,. 
The finite basis-size errors are still significant at this ba- 
sis size — at the CCSD(T) level, the binding energy is 
3.6 mE/j with the cc-pVQZ basis and 4.2 mE^ with the 
cc-pV5Z basis. Thus a more detailed study is necessary 
before a direct and definitive comparison can be made be- 
tween QMC and experiment (BE: 4.0 mE/j). The QMC 
result is consistent with the previous calculation using a 
plane- wave basis [9] , which should be at the infinite basis 
limit (aside from pseudopotential errors). 



2. Ionization potentials and electron affinities 



In Tables IV and V we present a summary of the cal- 
culated first ionization potential and electron affinity for 
first-row elements. For the ionization energy study, we 
used the cc-pVDZ and cc-pVTZ double- and triple-zeta 
basis sets [22] for all the elements. A few selected el- 
ements arc then studied, using the larger cc-pVQZ ba- 
sis sets. Also we looked at O with a cc-pV5Z basis 
set with 91 basis functions. For the electron affinity 
we used the aug-cc-pVDZ and aug-cc-pVTZ sets [46]. 
For comparison, we show the values calculated using HF, 
CCSD(T), and QMC, as well as the experimental data 
from Refs. [44, 45]. 



The agreement between QMC and the coupled cluster 
CCSD(T) results is in general very good. For unstable 
or meta-stable negative ions (N~ and Ne~), the Hamil- 
tonian expressed with a localized basis always yields an 
atomic-like ground state. In Tables IV and V, the total 
ground-state energies are not shown, but the mean dif- 
ference between the QMC and CCSD(T) values among 
all the atoms and ions is less than 3 mE^. The negative 
charged F ion, F~/aug-cc-pVTZ, has the largest discrep- 
ancy of ~ 7 mE,,. An MCSCF study [27] showed that the 
RHF single determinant gives a rather poor description 
of the ion. 
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TABLE IV: The first ionization potential of first-row ele- 
ments as calculated in the cc-pVDZ and cc-pVTZ basis sets. 
For selected elements, additional cc-pVQZ and cc-pV5Z ba- 
sis sets arc also used. Wo show the values obtained using 
HF, CCSD(T), and QMC, together with experimental results 
[44, 45]. QMC statistical errors axe in the last digit and are 
shown in parentheses. All energies are in eV. 



Atom 


HF 


CCSD(T) 


QMC 


Exp 


cc-pVDZ 










Li 


5.342 


5.345 


5.345(1) 


5.39 


Be 


8.079 


9.290 


9.296(2) 


9.32 


B 


8.038 


8.066 


7.942(5) 


8.30 


C 


10.803 


10.983 


11.025(5) 


11.26 


O 


11.965 


12.853 


12.808(2) 


13.62 


N 


1 ^ SQ^ 

LO.OijO 








F 


15.640 


16.710 


16.725(8) 


17.42 


Ne 


19.668 


20.893 


20.948(5) 


21.56 


cc-pVTZ 










Li 


5.342 


5.353 


5.353(1) 


5.39 












B 


8.038 


8.228 


8.08(1) 


8.30 


C 


10.798 


11.184 


11.237(8) 


11.26 


O 


12.011 


13.326 


13.256(7) 


13.62 


N 


13.892 


14.449 


14.574(4) 


14.53 


F 


15.654 


17.142 


17.154(7) 


17.42 


Ne 


19.673 


21.309 


21.37(1) 


21.56 


cc-pVQZ 










B 


8.041 


8.260 


8.19(5) 


8.30 


O 


12.018 


13.493 


13.44(1) 


13.62 


F 


15.649 


17.314 


17.33(2) 


17.42 


Ne 


19.661 


21.488 


21.57(2) 


21.56 


cc-pV5Z 










O 


12.022 


13.558 


13.50(1) 


13.62 



VI. H2O BOND STRETCHING 

We next examine the accuracy of our method in de- 
scribing bond stretching. The ability of a computa- 
tional method to deliver uniform accuracy as bonds are 
stretched /broken is obviously important in chemistry. It 
also provides an indicator for the potential accuracy of a 
method for solids, mimicking different levels of electron 
correlation. Our method uses a trial wave function in 
the constraint to deal with the phase problem. Almost 
all calculations to date, both with the plane-wave basis 
and in the present study, have used single Slater deter- 
minant trial wave functions. Since the quality of such 
wave functions decreases as bonds are stretched and cor- 
relation effects become more important, it is a significant 
challenge for the method to maintain its quality and ob- 
tain uniformly accurate results. 

Wc apply our method to stretched bonds in H2O, 
which has served as an excellent benchmark system 
[3, 32]. Symmetric 0-H bond length stretching is stud- 
ied. The bond angle between the two 0-H bonds is held 
fixed, while the 0-H bond length is increased from its 



TABLE V: Same as Table IV, but for the electron affinity and 
using the aug-cc-pVDZ and aug-cc-pVTZ basis sets [46]. 



Atom 




HF 


CCSD(T) 


QMC 


Exp 


aug-cc- 


-pVDZ 










B 




-0.300 


0.172 


0.125(5) 


0.28 


C 




0.468 


1.145 


1.205(6) 


1.263 







-0.523 


1.189 


1.25(1) 


1.46 


N 




-1.858 


-0.512 


-0.52(1) 


-0.07(2) 


F 




1.284 


3.228 


3.40(1) 


3.4 


Ne 




-7.724 


-7.321 


-7.32(1) 




aug-cc- 


-pVTZ 










B 




-0.305 


0.235 


0.23(4) 


0.28 


C 




0.453 


1.226 


1.29(3) 


1.263 







-0.565 


1.336 


1.40(5) 


1.46 


N 




-1.813 


-0.296 


-0.28(3) 


-0.07(2) 


F 




1.195 


3.317 


3.59(4) 


3.4 


Ne 




-6.479 


-6.108 


-6.30(7) 





equilibrium value. We considered three different basis 
sets: ST0-6G, cc-pVDZ [22], and DZ ANO [20]. For the 
small basis, we used the same geometry as in Table I. For 
the two larger basis sets, we used the same geometries as 
those of Ref. [32] and Ref. [3] for comparison with their 
FCI and DMRG energies, respectively. 

The results are shown in Table VI. In addition to 
QMC/UHF, i.e., using the variationally optimal HF so- 
lution as the trial wave function, we also carried out 
QMC calculations using the RHF solution, in order to 
further examine the effect of the trial wave function. As 
bonds are stretched, static correlation becomes increas- 
ingly important. The RHF solution becomes inceasingly 
unfavorable compared to the UHF solution, which has 
the correct behavior in the dissociation limit. As can be 
seen in Table VI, QMC with the RHF trial wave function 
(QMC/RHF) performs worse than QMC/UHF, which is 
consistent with the trend seen in Table III. Indeed, the 
QMC/RHF results deterioate as the bonds are stretched, 
reflecting the qualitatively incorrect nature of the RHF 
trial wave function at large bond lengths. 

In the range of bond lengths in Table VI, QMC with 
the UHF trial wave function gives roughly comparable 
accuracy as CCSD(T). For example, in the cc-pVDZ ba- 
sis at 1.5 i?e, QMC/UHF over-estimates the energy by 
2.7(1) mEh, while CCSD(T) over-estimates it by 1.6 mEh- 
At 2i?e, both QMC/UHF and CCSD(T) under-estimate 
the energy, by 3.0(2) and 3.8 mEh, respectively. In the 
DZ ANO' basis for 1.5 i?e, QMC/UHF is above the 
DMRG value by 2.1(5) mEh, while CCSD(T) is above 
by 1.6 mEh. 

Larger bond stretching of up to 8 Re is presented in 
Figure 4, using the cc-pVDZ valence double-zeta basis. 
QMC results are compared with exact FCI [32] and with 
coupled cluster results using RHF and UHF reference 
states (CCSD(T) and UCCSD(T), respectively). The in- 
set shows the errors of the various methods from the FCI 
numbers. CCSD(T) is excellent near equilibrium, but 
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TABLE VI: Basis-size and trial wave function dependence in stretched-bond calculations in H2O. Equilibrium geometries are 
the same as in Table. I. The all-electron total energy of H2O is shown using three different basis sets: minimal STO-6G, 
cc-pVDZ, and DZ ANO [21] with 7, 24, and 41 bases functions, respectively. The O-H bond lengths are stretched to 1.5 -Re 
and 2 Re. QMC energies obtained with both RHF and UHF trial wave functions are shown in the last two columns. Statistical 
errors in QMC are in the last digit, and are shown in parentheses. FCI/DMRG results are from the same sources as in Table I. 
All energies are in Hartrees. 

"bond length RHF UHF CCSD(T) FCI QMC/RHF QMC/UHF 
STO-6G 

1.5 Re -75.440 432 -75.502 069 -75.600 039 -75.576 8(3) -75.596 5(6) 

2_Re -75.141587 -75.464 541 -75.486 528 -75.355 7(3) -75.488 0(3) 

cc-pVDZ 

1.5 Re -75.802 387 -75.829 813 -76.070 717 -76.072 348 -76.068 8(3) -76.069 7(2) 

2 Re -75.587 711 -75.793 668 -75.955 485 -75.951665 -75.897 3(4) -75.954 6(2) 

DZ ANO 

1.5 Re -75.817 273 -75.852 670 -76.129 442 -76.131050 -76.126 5(7) -76.129 0(5) 

2 Re -75.602 850 -75.818 969 -76.009 395 -75.950(1) -76.011(1) 




I- 

1 2 3 4 5 6 7 

R/R 



FIG. 4: A comparison of bond stretching in H2 0/cc-pVDZ 
between the present QMC, coupled cluster methods, and 
FCI [32]. The main graph shows the calculated total energy 
(in Hartrees) as a function of the O-H bond length. The in- 
set shows the errors (in milli-Hartrees) of various methods 
with respect to FCI. Points are connected by straight lines 
for clarity. QMC statistical error bars are shown. 



is much worse at large bond lengths. (As mentioned, 
QMC with RHF trial wave functions would show similar 
behavior at large bond lengths.) UCCSD(T) performs 
much better for larger bond lengths, but is worse in the 
intermediate regime with errors larger than lOniEh. Our 
QMC results are in good agreement with the exact ener- 
gies, showing a maximum discrepancy of about 4(1) mEh, 
which occurs at 8Re. Together with the equilibrium re- 
sults, the method is seen to yield rather uniform accuracy 
across these bond lengths. This is encouraging, given 
that it is achieved with the same choice, namely the vari- 
ational UHF solution, as the trial wave function through- 
out the entire region. 



VII. DISCUSSION AND SUMMARY 

The present QMC method provides an approximate, 
but non-perturbative approach for many-electron calcu- 
lations. It can directly incorporate traditional electronic- 
structure machinery, such as high-quality basis sets, 
effective-core potentials, etc. The method obtains the 
many-body ground state by building a stochastic linear 
superposition of non-orthogonal independent-particle so- 
lutions. Thus, algorithmically it shares many of the in- 
gredients in more standard methods in quantum chem- 
istry and solid state physics. Its computational cost 
scales with the number of basis functions as to N'^. 

The accuracy of the present QMC method depends on 
the reference or trial wave function |4't), which in this 
study is chosen as the HF state. As in coupled cluster 
methods, the trial wave function is the starting point 
of the ground-state projection, but in QMC all excita- 
tions are implicitly included by the projection. Errors 
arise from the phaseless approximation in the projection, 
which uses |^'t) to approximate the phase of the contri- 
bution of a Slater determinant to the ground state, via 
the importance-sampling transformation. 

We have found, as shown in Tables HI and VI, that 
the most accurate energies are in general obtained us- 
ing the best variational single determinant |^t)- This 
is simply the HF solution when RHF and UHF are the 
same (e.g., in the H2O molecule at equilibrium), and is 
the UHF solution when the two differ. In the latter case, 
the method seems relatively insensitive, within the unre- 
stricted framework, to whether a Hartree-Fock or hybrid 
B3LYP Slater determinant is used. 

In contrast to the diffusion Monte Carlo (DMC) 
method, the present QMC must deal with finite basis-set 
errors and basis-size convergence. However, the ability 
to use any single-particle basis can also be advantageous. 
For example, in addition to the connection to quantum 
chemistry methods, our new method also makes possi- 
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ble QMC calculations for general model Hamiltonians, 
which are frequently used in the study of correlated elec- 
tron systems. Like the present AF QMC method, DMC is 
also approximate, using a trial wave function to impose 
the fixed-node approximation to control the sign prob- 
lem. The trial wave functions that have been used in our 
method are much simpler. Because of the complementary 
nature of the two different methods, direct comparisons 
are not always straightforward. Indications from the cur- 
rent calculations and earlier plane-wave studies are that 
the systematic accuracy of the phaseless approximation 
compare favorably with fixed- node DMC. 

In summary, we extended the recently introduced 
phaseless QMC method to handle any one-particle basis, 
and applied it to atoms and molecules using Gaussian ba- 
sis sets. Overall, our results at and near the equilibrium 
geometries are roughly comparable to those obtained us- 
ing CCSD(T), but are superior for bond stretching. Our 
preliminary results (to be published elsewhere) on bond- 
breaking in several diatomic molecules show similar uni- 
form accuracy. For a first application, these results are 



quite encouraging. There arc many possibilities for fur- 
ther improvement of the method. Currently, we are inves- 
tigating different forms of Hubbard-Stratonovich trans- 
formations, possibilities for better scaling, as well as ap- 
plications of the present method in transition metal ox- 
ides and other systems. 
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